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Abstract 

We compute the rapidity dependence of particle and transverse energy production in 
ultrarelativistic heavy ion collisions at various beam energies and atomic numbers using 
the perturbative QCD + saturation model. The distribution is a broad gaussian near 
y = but the rapid increase of particle production with the beam energy will via 
energy conservation strongly constrain the rapidity distribution at large y. 
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1 Introduction 



The Relativistic Heavy Ion Collider RHIC has already produced a lot of data at y/s 
up to 200 GeV and the planning for the ALICE/Large Hadron Collider experiment 
at y/s up to 5500 GeV is under way. The first data are on various large cross section 
phenomena like charged multiplicity at zero (pseudo)rapidity in nearly central colli- 
sions |1|, or at varying ^ impact parameter, or as a function of both rapidity 
and centrality ||, ||. Much theoretical effort, reviewed in [0, has been devoted both 
to predict the multiplicities before measurements and to draw conclusions from the 
completed measurements. With this information the reliability of predictions for the 
LHC energies, ^/s = 5500 GeV, is considerably enhanced. 

One of the models, reasonably successful in predicting the data at RHIC energies, is 
the saturation+pQCD model f^, based on microscopic 2^2 partonic processes. It is 
actually a member of a large class of models in which there is one dominant transverse 
momentum scale, Psat? determined by different versions of a saturation condition [1^, 
|TT], |T2|, |TB|. The purpose of this note is to apply the model to a study of the rapidity 
dependence of particle and transverse energy production. Related work in and 
T3[] , based on microscopic 2— >1 processes, will be compared with. 

The main observation is that, not surprisingly, energy conservation places strong 
constraints |T^ to the domain of validity of models based on independent 2— ^>2 subcol- 
lisions. There is no problem at zero rapidity: N{y = 0) and Exiy = 0) grow rapidly 
with ^/s, but still much more slowly that the total available energy ^/s. At larger 
rapidities such a rapid growth of N{y) cannot be sustained, the total energy carried 
by the produced particles would surpass the available total energy indicating that one 
has entered a new domain in which the independent scattering model no more is valid. 
We estimate that at RHIC the saturation+pQCD model could be valid up to ?/ = 1...2, 
and at LHC up to y = 3. ..4. Within this range the rapidity distribution is very flat, 
practically a wide gaussian with calculable curvature near y = 0. This is in marked 



contrast to the 2 — 1 model in |T1, ffH], where N(y) is sharply peaked at y = 0, in fact. 



has a discontinuous derivative there. 



2 Rapidity dependence in the saturation + pQCD 
model 

We shall first carry out an approximate analytic discussion emphasizing parametric 
dependences. One starts from the perturbatively determined leading-order two-jet 
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cross section (in standard notation): 



, 2^?° , =KY^ Xifi/A{Xi,p^) X2fj/A{X2,PT) ^ -JT ■ (1) 

dp^dyidy2 h dt 

To maintain the framework of @], we use the GRV94 leading-order parton distribu- 



tion functions [0, combined with the nuclear effects (shadowing) from the EKS98 



parametrization [|T^. An effective factor ^ 2 is to simulate the NLO contributions 



which are discussed in more detail in [^. After integration over y2 Eq. (^ gives a 
single-jet cross section behaving near the dominant saturation scale (which depends on 
^/s and y4, and also on y) parametrically like [ pi[ 

dyd^PT Pt 

where 5 ~ 0.5 (when no shadowing is included) and the ?/-dependence has been ap- 
proximated by a gaussian (the letter a is used for both the cross section and the width 
of the gaussian). A further integration over pt from some lower limit gives the hard 
cross section (see Figs. ^ and ^ below) 



dy pI 



Consider now central A+A collisions with the nuclear overlap function T^Aib = 0) = 
/{ttR^). The saturation condition then can be formulated as 

dN A? dapQCB^ . 2 

A"^ r:25a'l{pll y2/2a^po) 



7tR\ pI 

the solution of which gives po = Psat(^, y/s,y)- The saturation is here defined as the 
geometric saturation on the transverse plane of the two final state particles in the 
2— >-2 collision. The 2^1 models describe saturation as the separate saturation of 
the distribution functions of the two initial state particles. As long as there is one 
dominant tranverse momentum scale, the models lead to similar results, at least when 
final state particle production in nearly central collisions is considered. 

The principle underlying the saturation condition is illustrated by Fig. |I[ The per- 
turbative cross sections, of course, are valid only for large pt- However, one is only 
interested in the integral over pt, not the value of the distribution at small pt- The 
value of the full integral can also be reproduced by integrating the perturbative distri- 
bution from a lower limit given by the saturation condition. 
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Figure 1: The qualitative pT distribution of gluons produced in a heavy ion coUisions at 
^/s = 5500 GeV (sohd curve). Only the integral over px has physical meaning; the satura- 
tion+pQCD model produces the same integral by integrating over the dashed curve computed 
using pQCD but with the lower limit determined by the saturation condition (|5|). 



To approximately solve Eq. (H) for pq = psat one may neglect the po-dependence of 
a and approximate <y'^{po) ~ VPo^i ^ ~ 0.44. This numerical value actually contains 



also a contribution from the p-p- integration [^. Then, 

~ 0.208 GeVA°-i28y^o.i9i^_,2/(4(2+^)^.)^ 
dN/dy ~ A(6+^«/(6+30^2V{2+o^-,V(2(2+0^^) (7) 

where the numerical values computed in with shadowing are also given. 

The scaling exponents of A and ^/s are the ones from p, |21| ; the rapidity dependence 



is new. Note how the saturation condition leads to ~ A"-, a ~ 1 [|Il|], though in 
the independent collision limit ~ A^^^. Similarly, the saturation condition has a 
significant effect on the rapidity distribution of saturated particle production: it is 
wider than that of elementary subprocesses by a factor y/2 + ^ ^ 1.6. 

In the numerical computation, one first evaluates dN/dy{po) = TAA{0)d(TpQCD / dy{po) , 
The result is shown in Fig. ^together with a gaussian fit ~ exp(— ?/^/(2(T^)). The dis- 
tributions get narrower when pq increases: at y/s = 5500 GeV cr(po) ~ 6.6 — po- 
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Figure 2: Plots of dN/dy = TAA{0)dapQCD /dy{po) for ^ = 200 GeV (from top to bottom 
po = 1.0,1.2,1.5,1.7,2.0 GeV) and 5500 GeV (from top to bottom po = 1.5,1.7,2.0,2.2,2.5 
GeV). The dashed hnes show a gaussian fit ~ exp(— y^/(2cr^)). 



Using the numerically computed dN/dy{po) the saturation condition operates as 
shown in Fig. ^ and leads to NaaIu) also in Fig. ^. In agreement with the analytic 
argument in Eq. (Rl) a very broad distribution near y = is obtained: 



cr 



5.8 
5.9 
6.8 



's = 130 GeV, 
= 200 GeV, 
^ = 5500 GeV. 



(9) 
(10) 

(11) 



The form of the rapidity dependence of the saturation scale proposed in ||14| , p!5| is 
markedly different. One first notes that the conjectured saturation scale as determined 



from deep inelastic scattering [22| can be parametrised as 



Q^(x) = Ql i^) 



X J 



(12) 



where Qo ~ 1 GeV, ~ 3 ■ 10 ^ and A ~ 0.3. To relate the DIS Bjorken variable 
X to the A+A gluon production variables, one further writes that x ~ e 



y-VB 



Then 

Q'siy) ~ Aj/i^e"'^'^'. Including only the region in which both initial gluons in the 2^1 
process are saturated, one obtains 



^ - Qliy) 



(13) 



Thus the rapidity distribution will have a sharp peak at y = 0, as shown in the right 
panel of Fig. 0. In Ref. |1^ it is assummed that the pseudorapidity distribution of 
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Figure 3: The left panel shows how the saturation condition (^) operates: the dashed line 
shows its RHS PqR\, the solid lines show the LHS T^yi(0)(i(TpQCD(P0) y)/dy at ^/s = 5500 GeV 
for various y as a function of po- The point of intersection gives both the value of psat(y) and 
of dN / dy{psat{y)) , the value of which is plotted in panel at right for central Pb-Pb collisions 
at -^/i = 5500 GeV and Au-Au at 200 GeV. The dotted and dashed curves show the prediction 
of for A = 0.25, 0.3 respectively. 



produced gluons (virtual, ^ Qg ■ 1 GeV, pr ^ Qs) directly give the pseudorapidity 
distributions of pions in the final state, apart from a multiplicative normalization factor 
c. Thus dNch/dr] ^ c ■ dNg/dr], where c contains factors from the actual gluon produc- 
tion, such as the gluon liberation constant Cg = 0{1) and possible factors related to 
the normalization and approximations made for the unintegrated gluon densities. The 
constant c also includes factors related to hadronization, decays and particle content 
in the final state. Such factors are c-^/g which indicates how many pions come from 
each virtual gluon, and a factor Cch accounting for the conversion of total to charged 
multiplicity. It should be noted that when plotting the rapidity distributions of [1^1 for 
the produced gluons in Figs. 3 (right panel) and in Fig. 5 (left panel), we have fixed the 
normalization constant c = Cch ■ c^^/g ■ Cg based on the measured dNch/dr] and used 
Cch ~ 2/3. As our primary goal here is to compare the shapes of the rapidity spectra 
of the two approaches, we have not tried to undo the effect of the factor c-j^/g- As the 
virtuality of the produced gluons in |14| is taken to be larger than rrip, we expect that 
the actual number of initially produced gluons in the saturation model of |T^ can be 
down by c^r/g ~ 2 relative to what is shown in Figs. 3 and 5. 

One may observe that in the final state saturation calculation one does not at all 
need the distribution functions within the saturation region as defined by Eq. ([1^). To 
check the effect of saturated distribution functions on the 2^2 model one can perform 
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the following computation. Including only the dominant gluons one has 



da 



dyd'^pj 



9ai 



(14) 



where Xi = Pt/ ■ (e^ + e^^ 
distribution functions (see 



X2 = Pr/y/s- (e ^ + e ^2). A simple model for saturated 
) would be 



Jb 



(15) 



S ( 

xg{x,PT) = C— \ pt-Q 

with 5*^ ~ Tri?^, C = constant, and where implicitly Qq ~ A^^^. Integration over j/2 and 
Pt then gives a ^/-distribution which is broad and very similar to those in Figs. ^ and ^. 
We thus conclude that the broad y-distribution is not due to the use of distributions 
in the DGLAP region, it is a property of the 2^2 model as such. 

The a/s and y dependences of the result (|T3D are sensitive to the value of the param- 
eter A = 0.25 . . . 0.30, especially when one approaches LHC energies. Furthermore this 
A interval is based on the analysis in Ref. |^ which neglects the scale evolution of the 
gluon distribution in the non-saturated region on which the saturated gluon distribu- 
tion has to be matched in the vicinity of the 'critical line' Ql{x) defining the transition 
region from saturated to non-saturated gluon densities. An interesting recent analysis 
in 1^ suggests that one might be led to consider values as large as A ~ 0.5. This 
is similar to the behaviour of the gluon distribution at the final state saturation scale 
2l. 



3 Energy conservation 

The results (^)-(H) lead to a rather striking powerlike growth with energy, especially 
for the transverse energy ~ Psatiy)N{y). Let us thus estimate when the total amount 
of energy within some rapidity range — ^up < y < 1/up is less than some fraction, 40%, 
say, of the total energy Ay/s ~ Ae^^ available (GeV units are used), neglecting first 
the weak rapidity dependence. A simple analytic estimate of this is (cat, Cp are the 
constants in (|D,(P)): 



PsAy = 0) coshy ^ — 
y=o dy 



ry^p , dN 
-yup dy 

-3S/{2+0 



Psat(0)e^"P (16) 

y=0 



ry up 

E{\y\ < Vup) ~ / dy 
J-yup 

giving 

2/,p + log(c^c,A(i-«)/(6+30) < ^^(1 _ J^) + iogo.4. (17) 
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The two logs are constants of order one, but the qualitatively important factor is 
— 35/(2 + ^)) ^ — 0.62); it expresses the fact that with increasing beam 
energy particle production in the saturation + pQCD model will exhaust the total 
available energy within a fixed fraction of the beam rapidity ub- Beyond that rapidity 
this independent subcoUision model cannot be valid and the physical mechanism has 
to change completely. The same effect is shown more quantitatively in Fig. ^ where 
we have computed the total produced energy as 

/■yup f dN 
E{y^p) = dy dpr-, — ^--pTCOsh?/, (18) 

The rapidity distributions of produced partons are obtained from Eq. (|l]), as discussed 
in detail in Ref. |2^. Keeping the saturation scale fixed at Psat{y = 0) gives the 
dotted lines, not qualitatively different from the solid lines obtained with the rapidity 
dependent saturation. 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 



Figure 4: The total energy produced in the saturation + pQCD model up to \y\ < yup 
relative to the cms-energy in central A A collisions. The cms-energies and nuclei are denoted 
in the figure. "Initial" refers to time after formation before any collective expansion. The 
solid curves are computed from Eq. (|^) with Psat(y)i and the dotted curves with a fixed 
Psax{y = 0) as the lower limit of pr-integration. 

The saturation-|-pQCD model thus predicts that the ?/-distribution Niy) is very flat, 
gaussianlike around ?/ = 0, with height increasing as given by Eq. (|^) but extending 
only up to some maximum value of y imposed by energy conservation. We emphasize 
that the saturation condition itself leads to a broadening of the rapidity distribution. 
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Allowing 40% of the total energy leads to \y\<2 at RHIC (with ys = 4.93 (5.36) for 
= 130 (200) GeV) and \y\<3.5 at LHC (with ys = 8.67). Beyond these values 
correlations between subcollisions must start to play an increasingly important role. 



4 Hydrodynamic evolution 

The previous considerations apply at the initial time Ti{y) ^ l/psa.t{y)- Within the 
domain of validity, \y\ < yup, there is little variation in Ti{y). Hydrodynamic evolution 
in the case of y-dependent initial conditions of the type obtained here was numerically 
studied in p5|. There is always significant flow of Ej- from small to large rapidities, 
due to pdV work p, ^] , but in this y-dependent case there is also some flow of entropy 
from small to large rapidities. Typically, the entropy density at y = decreases by 
about 10%, as shown in . 



We thus conclude that the time evolution of y-distribution is one more detail to be 
added to the list of factors affecting the relation between initial and final multiplicity 



26|| : initial and final particle-to-entropy ratios, charged-to-neutral ratio, transverse 
expansion, equation of state, decouphng effects and resonance decays. However, in the 
central region the hydrodynamic evolution affects only the overall magnitude of the 
rapidity distribution while the effect on the shape is small. 



5 Comparison with RHIC data 

We note first that the RHIC data 0, ^ gives the pseudorapidity distribution dNch/dr] 
while our computations so far are for dN/ dy at the moment of formation of the partonic 
system. The conversion of to the observed Net is simple (see end of previous section). 



but the relation between y- and ?7-distributions is more subtle ^6 
The rapidity distribution of particles is the sum 

dN r , ^ dN, 



f Ctl\j , , 

(19) 



dy J dprdy 

where the index i runs through all particles included in the spectrum, similarly for 
pseudorapidity. We have shown explicitly the transverse momentum integration to 
emphasize that the connection between the rapidity and pseudorapidity depends on 
the mass and the transverse momentum of particles: y = SiTsmh.{{pT/mTi) sinhr)). 
Similarly the Jacobian 

Ji{r],pT,mi) = ^ = (20) 

in the transformation between the spectra 

dN _ f dN f ^ dNi{pT,y) 
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depends on the mass and the transverse momentum. In the following we will assume 
that pions dominate and consider equations for single particle only. 

The pseudorapidity data |^ is characterized by a dip in a central plateau of total 
extent of 3... 4 in 77. Starting from a flat rapidity distribution, this dip results from the 
Jacobian in the transformation to pseudorapidity. The depth of the dip is quite sensitive 
to particle masses and shapes of the transverse momentum and rapidity distributions, 
variations by a factor two in the depth of the dip can easily take place. We also empha- 
size the difference between performing the transformation using Eq. ( pT]) with px and y 
distributions or simply approximating the Jacobian by dN/dr] = J{ri, {pT),m)dN/dy. 
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Figure 5: (Left) Initial rapidity distributions dN/dy for the saturation+pQCD model (solid 
curve, EKRT), for the double- Woods-Saxon (dWS) parametrization (|22|) (dotted curve) and 
for the prediction of |0 (dashed curve) (see text). (Right) The final state pseudorapidity 
distributions dN/dr] (6% central): data 0, the saturation + pQCD model (thick solid lines, 
EKRT) and the parametrisation ( p2|) (see text). 



In Fig. I^a the results from calculations of the initial dN/dy are shown for nearly 
central Au+Au collisions at ^/s = 130 GeV. The solid curve (EKRT) is the prediction 
with y-dependent saturation of produced partons, computed as in Figs. ||. A 6% 
centrality cut corresponds to an effective nucleus A = 178, as discussed in Within 
the dashed part of the curve more than 40% of total energy has been consumed (Fig. H); 
the curves stop when all the energy is consumed. The dotted curve shows a double- 
Woods-Saxon (dWS) parametrization 



(1 + e-y°^'^f 



I _|_ e(-y-yo)/rf)(i -\- g{y-yo)/d^ ' 



(22) 
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normalised to the EKRT result at y = 0. Here yo = 3.3 gives the width of the 
distribution and d = 0.65 the steepness at this cutoff. This double Woods-Saxon form 
for dN/dy is an arbitrary parametrization designed to reproduce the pseudorapidity 
data but it has the two elements we want to emphasize, the flatness of dN/ dy at small 
rapidities and the rapid decrease at large values of y required by energy conservation. 



Finally, the dashed curve is the initially produced rapidity distribution from [|14[. The 
normalization of this curve was discussed in the end of Sec. 2. 

Fig. shows the data with dNch/dr] curves calculated from dN/dy assuming 
that all particles are pions with an exponential transverse momentum distribution, 
dN/dpldy ~ dN/dy ■ exp(-2pr/(pr)), with {pr) = 0.4 GeV or 0.5 GeV. The thick 
solid lines (EKRT) are the saturation model prediction, computed using Eq. (^1]). 
The dashed parts of the lines correspond to those in the left panel and indicate again 
where energy conservation is expected to suppress the distributions. The thin solid 
lines (dWS) are from the parametrization (^) with yo = 3.3 and d = 0.65, which 
reproduce the PHOBOS data. In the framework of P], the initial and final state 
rapidity distributions are connected through entropy conservation by dNch/dy = 2/3 ■ 
3.6/4.0 ■ dN/dy{psa,t{y))- Here, in order to study the shapes of the obtained spectra, all 
the curves are normalized to the same point at r/ = 0. Compared with the normalization 
in §, an additional factor 0.95 (0.92) for (pr) = 0.4 (0.5) GeV is applied to the EKRT 
curves in Fig. ^d. The relation between the initial and final multiplicities is studied in 



more detail in . 

We do not show here the pseudorapidity distribution from the dashed curve in the 
left panel; in WM the y ^ t] transformation was carried out simply using a Jacobian 
J{ri,pT,mg), calculated for virtual gluons with ^ Qg ■ 1 GeV and pr ~ Qs, as a 
multiplicative factor between dN/dy and dN/dr]. With Qsiy = 0) = -\/2 GeV the 
suppression from the y- to t]- distribution at y = rj = is pt/ttlj' ^ 0.76, which turns 
the peak in the y distribution into a dip in the //-distribution. The results of are 
thus very close to the data in Fig. ^d. 



6 Conclusions 

We have discussed the rapidity distribution dN/dy of particles produced in very high 
energy central A + A collisions in the saturation+pQCD model. The distribution 
around y = will be a broad gaussian with a calculable width and with the value 
dN/dy{y = 0) increasing rapidly, powerlike. However, the model can be applied only 
until some value y^^p (~ 1...2 at RHIC, 3. ..4 at LHC), beyond which some new type of 
fragmentation region dynamics, taking into account correlations between subcoUisions, 
must enter. At asymptotic energies the saturation region is parametrically favoured by 
powers of I /as and this will damp the distribution at large y. 

The overall magnitude of dN/dr] at = predicted by the saturation + pQCD 
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model agreed very well with RHIC experiments for central collisions. The centrality 
dependence was also reproduced on the 10% level within A^part^lOO. A more reliable 
calculation on the ?/-dependence of initial gluon production can be made only after a 
better understanding of the correlations among subprocesses in the fragmentation and 
near fragmentation regions has been achieved. At present one may note that within 
the domain of validity of the model, before energy conservation makes the assumption 
of independent subcollisions invalid, there is agreement with data. The width of the 
relevant y-range at RHIC is narrow, but becomes appreciable at LHC energies. It 
will be very interesting to see how a rapid increase at y = and dynamics in the 
fragmentation regions will be related. 
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